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Abstract 

We present an algorithm (CaSPA) which accounts for the effects of 
periodic boundary conditions in the calculation of size of percolating ag- 
gregated clusters. The algorithm calculates the gyration tensor, allowing 
for a mixture of infinite (macroscale) and finite (microscale) principle mo- 
ments. Equilibration of a triblock copolymer system from a disordered 
initial configuration to a hexagonal phase is examined using the algorithm. 
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1 Introduction 

Structures obtained from surfactant solutions or block copolymers and their 
properties are of interest in many fields jl]. The processes that occur during 
these transitions can be relevant in polymer processing for the fabrication of 
nanostructured materials [2], especially in cases where the reversibility between 
structures is of interest. In biophysics, the lamellar to reverse hexagonal phase 
transition is considered as the first step in understanding membrane fusion [3l 

Order-disorder and order-order transitions have been studied experimentally, 
theoretically and using computer simulations. Much of the theoretical and simu- 
lation work has focused on the identification of different ordered structures that 
can be obtained when changing the architecture of the amphiphilic molecule and 
the system conditions [6]. Nevertheless, the dynamic processes that describe the 
order-disorder or order-order transitions are also important. Processes such as 
micelle formation and stabiHzation, sphere to rod transitions, bilayer breakdown 
and structural changes during the formation of nanostructured materials have 
attracted much attention [7l [8] . 

Groot and Madden used dissipative particle dynamics (DPD) simulations to 
describe the formation of an hexagonal phase from a disordered phase, where 
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an unstable gyroid phase appears as an intermediate using DPD simulations [9]. 
More recently, Soto-Figueroa et al. described the dynamics of different order- 
order transitions in polystyrene-polyisoprene diblock copolymers [lO] . The tran- 
sition between hexagonally packed cylinders to an array of body centred cube 
spheres is a result of undulations in the cylinders that eventually break into 
ellipsoids to latter form spheres. They also describe the transition from a bi- 
continuous structure to a lamellar phase going through an intermediate phase 
containing infinite cylinders, not observed by Groot and Madden, before a lamel- 
lar phase is obtained. Dynamics of the formation of ordered phases has also been 
reported for a variety of surfactant architectures [TT], but most of the results are 
limited to a collection of snapshots at different times during the simulation. In 
some cases, an order parameter is defined and used to determine the evolution 
of the observed phases, which requires the calculation of the structure factor 



For any simulation approach which seeks to model mesoscale aggregates 
from the microscale, the aggregates will appear infinite, that is, they will per- 
colate, spanning the periodic boundary conditions (PBCs) of the simulation. 
The dimensionality of the aggregate (whether one dimensional for cylinders, 
two dimensional for lamellae, or three dimensional for network structures) will 
set the dimensionality of this percolation. Algorithms exist to identify aggre- 
gates within a simulation configuration (principally, the Hoshen-Kopelman (H- 
K) algorithm [l3]), however, once identified, the aggregate must be properly 
characterised. 

The standard approach to characterisation of the size and shape of an ag- 
gregate is diagonalization of its gyration tensor S. If the positions of the TV 
particles in the aggregate relative to its center of mass are given by {r^}, the 
gyration tensor is given by: 



where rf denotes the a'th component of the vector r, and Sap indicates the 
(a, /3)'th element of the tensor S. The eigenvectors of this tensor (the "principal 
axes") give the orientation of the aggregate, and the associated eigenvalues (the 
"principle moments") give the length scales of the aggregate along these vectors. 

A naive implementation of Eqn. [1] for a percolating aggregate, using only 
the coordinates from the microscale aggregate identified within the simulation 
boundary conditions, will produce incorrect results. First, such an implementa- 
tion will not identify the correct principle moments; since the full, macroscale 
aggregate is percolating, at least one principle moment is infinite by defini- 
tion. Second, there is some subtlety as to which periodic images of a particle 
to include in the calculation, such that only particles within a single aggregate 
should be included, and particles within images of the aggregate must be ex- 
cluded. Finally, by not expHcitly dealing with the percolation of the aggregate, 
pathological configurations may result in an incorrectly oriented gyration ten- 
sor. If the principle axes are incorrectly oriented, then equivalently the principle 
moments will be wrong. This is illustrated in[TJ 

In this work, we present an algorithm (Calculation of Size of Percolating Ag- 
gregates, or CaSPA) to deal with these issues. The algorithm is demonstrated 
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Figure 1: An aggregate (filled area) repeated across periodic boundary con- 
ditions. The micro-scale aggregate identified from the structure is shown by 
the black filled area, and the gyration tensor for this aggregate is schematically 
illustrated by the full arrows. To calculate the gyration tensor for the true per- 
colating aggregate, only particles in the black and dark shaded areas should be 
included. The gyration tensor for the percolating aggregate is schematically il- 
lustrated by the dotted arrows. The two gyration tensors are not aligned. Note 
that the percolating aggregate is infinite in the horizontal direction, whereas 
the longest principle component of the microscale aggregate lies closest to the 
vertical direction. 

through application to aggregates formed in a triblock copolymer system, sim- 
ulated using dissipative particle dynamics (DPD) [T4 t \TE \ flB t 117). 

2 Methods 

2.1 CaSPA algorithm 

Consider an aggregate composed of N particles {r.i}, with gyration tensor S. 
The aggregate exists within PBCs, and is a percolating cluster, that is, particles 
within the aggregate are connected to particles in certain periodic images of the 
aggregate. We denote the set of translation vectors between the contacting 
images {p^}; these can be reduced to a combination of Up linearly independent 
vectors {pi}. The aggregate is therefore percolating in Up dimensions. We seek 
the gyration tensor 8™'"^''° for the full macroscale aggregate, made up of the 
combination of individual, contacting aggregates across the periodic boundaries. 
If we include Ui image aggregates along each of the p vectors, this gyration tensor 
will be given by: 

/ («,-l)/2 K^l)/2 \ 

« \fci = -(„.-l)/2 fc„p=-K-l)/2/ 
N / rip \ / Up \ 

E ^? + Efc.p" ^f + Evp^' 

Multiplying out, through the symmetry of the sum limits, and considering 
that since the origin is at the center of mass X^iPf = 0, we find: 
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where the tensor Q is given by: 



5X'^° = 5a/3 + ^(n.-l)K + l)Q„^ (3) 



u = Y.p"p' (4) 



4=1 

Note the similarity between the form of Q and the definition of the gyration 
tensor (Eq. [T|). The tensor Q can be considered as the normahzed gyration 
tensor for a coarse-grained representation of the macroscale aggregate, with one 
point mass per image. We seek to diagonaHse 8'"°'^''°, giving three eigenvec- 
tors (the principle axes) , and three corresponding eigenvalues (the principle 
components) A^, < A2 < A|. 

As rii tends to infinity (the bulk limit), the gyration tensor 8™°'^''° becomes 
dominated by the matrix Q. This determines the orientation of the macroscale 
aggregate in space. However, provided the aggregate does not percolate in all 
dimensions, one or more of the eigen?;a/ttes of Q will be zero, hence Q is singular. 
By performing singular value decomposition [18j on Q, we can identify the range 
U of Q, the set of Up eigenvectors with non-zero eigenvalues, and the nullspace 
V of Q, a set of orthogonal eigenvectors for which the eigenvalues are zero. The 
macroscale aggregate will have infinite principle components in the range (where 
Q will dominate), and finite principle components in the nullspace (where only 
S contributes). 

For tensors Q with a single vector Vi in the nullspace, the orientation of 
the finite principle component of the macroscale aggregate is given by Ai = Vi 
and the value of the finite principal component is given by the projection of S 
along this vector, A^ — Vf SVi. The remaining two principle components are 
infinite, A2 = A3 = 00, with the corresponding principle axes lying along the 
range, A2 = U2 and A3 = U3. 

For a tensor Q with two vectors (Vi and V2) in the nullspace, the orientation 
of the infinite principle component A3 = cxd of the macroscale aggregate is given 
by A3 = U3. The vectors describing the finite principal components will be 
linear combinations of Vi and V2, which will be a pair of vectors which are 
orthogonal to the range U3, but are otherwise arbitrary. To find the finite 
principle axes, we must project the matrix S into the plane given by V. This 

\ VT ] r ^ 

results in a two dimensional gyration tensor given by _ .7- S [ Vi V2 J , 

with two-dimensional eigenvectors A'^ and A2, and corresponding eigenvalues 
giving the principle components, A^ and Aj. The orientations of the principle 
components are given by A^ = A^ [Vi V2], that is, the eigenvectors AJ- give the 
appropriate linear combinations of the nullspace vectors to give the principle 
axes. 

The two remaining cases are tensors Q with no nullspace, representing a 
percolating macroscale aggregate which is infinite in all directions (A^ = A2 = 
A3 = 00), and tensors Q with a three-dimensional nullspace, representing ag- 
gregates with no self contacts, which can be treated entirely from the aggregate 
gyration tensor S. 

From this, the algorithm to calculate the principle components and axes of 
the macroscale aggregate must perform the following steps: 
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1. Identify the aggregate coordinates {r^}, and calculate the aggregate gyra- 
tion tensor S. 

2. Identify the set of self-contact vectors, {p^}, and reduce it to a linearly 
dependent set, {pi}. 

3. Calculate the tensor Q ({Pi}), and find its range U and nullspace V. 

4. Calculate the projection of the gyration tensor S onto the nullspace V. 

5. Return the macroscale aggregate principle axes (the range U and unit 
vector projections of the nullspace V onto the gyration tensor S) and 
principle components (infinity for axes corresponding to the range, and 
the length of the vector projections of the nullspace V onto the gyration 
tensor S for the remaining components). 

We now describe how such an algorithm may be implemented. 

2.1.1 Identification of the aggregate coordinates 

The H-K algorithm is suitable for initial identification of an aggregate from 
a particle configuration. However, to deal with an aggregate which crosses 
PBC's, an extra step is necessary. In simulation coordinates, aggregates will 
often consist of a number of disjunct segments, connected across the PBC's 
(see [2]). The first task of the algorithm, once the initial iV-particle aggregate 
has been identified, is to "stitch" these disjunct segments into a single, fully 
connected object. The minimum image convention cannot be used for this, 
since the positions of particles in images of the aggregate may become mixed 
with positions of particles in the aggregate during calculation of the gyration 
tensor, giving incorrect results. 

To "stitch" an aggregate together, the coordinates {r^} of the particles in 
the aggregate are passed to a version of the cluster identification algorithm 
which does not recognise the PRCs. This identifies the disjoint segments, 
and returns a list {k}, labelling each of the N particles in {r^} according to 
which subcluster that particle is a member of. The algorithm then loops over 
all pairs of particles, until it finds a pair such that k ^ Ij, and which 

contact across PBCs. Having found such a pair, all particles with label U are 
relabelled with Ij and translated such that particles i and j are bonded without 
PBCs (merging disjunct clusters k and Ij into a single, connected cluster), and 
Hd is decremented. This is then repeated until — 1, that is, there exists 
only a single connected cluster. Note that the resulting set of coordinates will 
now extend outside of the original PBCs (see[2|). A fiowchart for the "stitching" 
algorithm is shown in Fig. [3]. Given this single connected cluster (the microscale 
aggregate), the gyration tensor S can be calculated as Eqn. [TJ 

2.1.2 Identification of the Self-Contact Vectors 

Once the aggregate coordinates have been identified, the self-contact vectors can 
be found. The set {p^} is easy to identify - loop across every pair of particles 
in the aggregate, check if they are contacting across the PBCs, and if they are, 
store the vector connecting the periodic images (given, in the terminology of [31 
by Shift{Ar)). This set may, however, be very large, and will in general contain 
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Figure 2: Illustration of the "stitching" process. Shaded circles represent parti- 
cles within the PBC's which form two disjunct segments of a single aggregate. 
Once self-contact vectors have been identified, one of these segments is trans- 
lated outside of the PBC's (open circles) such that the aggregate forms a single, 
connected cluster. 
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Figure 3: The "stitching" algorithm. Inset: the "shift" algorithm, which finds 

the vector between periodic images. The function rint rounds its argument to 
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a large quantity of redundant data. We seek instead the linearly independent 
set {pi}, which in three dimensions may consist of a maximum of three vectors. 

To find these vectors, the algorithm must loop across every pair of particles 
in the aggregate, checking for contact across PBCs, until a vector connecting 
periodic images is found. This is stored in pi. The loop then continues, until 
a second candidate vector p' is found. The algorithm must check that this new 
vector is linearly independent of pi; the condition for this is that |pi x p'| ^ 0. 
If p' is linearly independent of pi, then the algorithm sets P2 = p'. Otherwise, 
the algorithm continues to check candidate vectors p' for Hnear dependence 
against pi. 

If two linearly dependent vectors have been found, then the loop once again 
continues, but now candidate vectors must be checked for linear independence 
against pi and P2. The condition for linear independence is now that det ([pi P2 p']) ^ 
0. If p' is Hnearly independent of pi, the algorithm sets pa = p', and terminates; 
the aggregate is percolating in three dimensions, hence all principle components 
are infinite. 

Once all pairs have been checked, the algorithm is left with rip linearly inde- 
pendent self-contact vectors {pi}. A fiowchart for the self-contact identification 
algorithm is shown in|4l From the vectors {pi}, the tensor Q may be calculated 
as Eqn. [H This tensor is then diagonaHsed. There should be non-zero eigen- 
values, corresponding to the eigenvectors {Ui} in the range of Q, and (3 — rip) 
zero eigenvalues, corresponding to the eigenvectors {V^} in the nullspace of Q. 

2.1.3 Calculation of Gyration Tensor of the Infinitely Repeated Ag- 
gregate 

Finally, the algorithm has the aggregate gyration tensor S, the dimensionality 
of percolation, np, and the range and nullspace of percolation, {Ui} and {V^}. 
The algorithm will already have terminated if the aggregate is percolating in 
three dimensions, returning three infinite principle components. Otherwise, the 
principle components and axes may be calculated by projection of S onto the 
nullspace V, as detailed above. A full fiowchart of CaSPA is shown in[5l 

2.1.4 Optimization 

The algorithm as described here has been separated into parts for conceptual 
clarity, and as such is not optimized for speed. It is anticipated, however, that 
the generation of uncorrelated configurations by simulation for analysis is likely 
to be significantly slower than the analysis of these configurations. 

The process of checking all pairs in the system means that, as presented, the 
time needed scales as iV^. If the definition of contact allows, the algorithm can 
be reduced to scaling with N via construction and use of a cell list. A second ob- 
vious improvement would be to perform the identification of self-contact vectors 
{pi} at the same time as the "stitching" algorithm, by checking for boundary- 
crossing contacts between particles in the same aggregate segment. While this 
will not change the scaHng form, it should significantly reduce the constant of 
proportionality. 
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Figure 5: The CaSPA algorithm; "CoM frame" denotes the frame of reference 
where the aggregate center of mass is the origin. 
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2.2 Test System 



To demonstrate the algorithm, we have simulated equilibration of a triblock 
copolymer solution from an initial, random configuration. The simulation has 
been performed using dissipative particle dynamics (DPD). This is a mesoscale 
simulation methodology, where "particles" represent loosely-defined "fiuid ele- 
ments", and move via a combination of dissipative, random and conservative 
interparticle forces. 

Particles in the system interact via the "standard" DPD force field [16], given 



where r.^j is the separation between particles indices i and j, is a notional 
particle diameter (setting the length scale of the system) , and aij a parameter 
which depends upon the types of particles indices i and j, setting the strength 
of interaction between these types. Particles which are bonded along a polymer 
chain also interact via a spring force, spring constant C = 'iksT |17j. 

The simulated system exists within a cubic box of size 15rc, with total 
particle number density p — 3, and is made up of A and B type monomers 
and solvent particles (type S). Interaction parameters for like particle types 
are axx = 25, and for unlike types are as a — clab — 50 and asB — 25; as 
such, A-type monomers are "hydrophobic". All A and B type monomers exist 
as part of A3B10A3 triblock copolymers. The system has 4061 solvent particles, 
and 379 copolymers, giving a copolymer volume fraction of 0.6. The simulation 
reduced timestep is At — 0.01, and the simulation is run for 4 x 10^ timesteps. 
The initial configuration is generated by placing solvent particles and one end 
monomer of each copolymer randomly in the simulation box. Copolymers are 
grown by placing the center of the next monomer in each polymer at a random 
position on a sphere of radius Tc centered on the previous monomer. The energy 
of this random configuration is then minimised by steepest descent to give the 
starting configuration for the simulation. 



Time series data of the configurational energy (decomposed into pair and spring 
interactions) from the simulation are shown in [6l From the point of view of en- 
ergy, the random initial system configuration appears to refiect the equilibrium 
state well, and the system does not appear to pass through any significant en- 
ergetic relaxation. 

"Stitched" snapshot configurations from the simulation (seelTj) indicate that 
the system undergoes mesophase separation, with the "hydrophobic" ^-type 
monomers forming elongated aggregates, demixed from the solvent and B-type 
monomers. As such, the CaSPA algorithm has been used to study these ag- 
gregates. Pairs of A- type monomers are considered to be "connected" if they 
overlap (that is, if < Tc). Data for clusters smaller than 3 monomers are 
discarded. It should be noted that the results we present are insensitive to 
the definition of connectivity, being qualitatively similar across a range of def- 
initions of connectivity from r < O.Src to r < 1.2rc. This insensitivity can 
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Figure 6: Configurational energy for the triblock copolymer simulation, de- 
composed into pair interactions (C/c) and spring force interactions (Us), as a 
function of simulation time t. No significant energetic relaxation is observed. 




Figure 7: "Stitched" snapshot conformations from the triblock copolymer sim- 
ulation. Cylinders join A-type monomers separated by less than Vc- Lines 
indicate the simulation box. a) Coexistence between large, self-contacting ag- 
gregate (lower right) and elongated micelle at t = 2At. b) Coexistence between 
large, self-contacting aggregate (center), spherical micelle (top) and elongated 
micelle (bottom) at t — lOOAi. c) Single, elongated micelle at t = 200Ai. d) 
Large, self-contacting aggregate at t — 380Ai. 
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Figure 8: Radial distribution function gAA{r) between hydrophobic monomers. 



Non-percolating Percolating in ID 

20 h' ' ' ' \ ' ' ' ' ' k ' i ' ' ' 'J I ' \ ' \ ' \ ^ 




^ ' 



' ' ' ' - 

10000 20000 30000 40000 10000 20000 30000 40000 

t/At t/At 



2 

n 



Figure 9: Principal components of the gyration tensor Ai < A2 < A3 for clusters 
of ^-type monomers as a function of simulation time t. The left hand column 
shows principal components for non-percolating clusters. The right hand column 
shows principal components for clusters which are percolating in one dimension; 
note that, for these clusters, A3 = cx) by definition. Multiple points at a given 
time step denote multiple clusters of that type at a given time step. 



be explained by considering the radial distribution function between A-type 
monomers gAA{r) (see [8]). The nearest-neighbour peak can be seen to be at a 
maximum at at r = 0.825rc, and can further be seen to be distinct from the 
broader peak indicating aggregate structure up to r « 1.2rc. As long as the 
cutoff for connectivity lies between these two values, the algorithm will probe 
nearest neighbour structure in an effective fashion. 

Time series data from CaSPA for the principal gyration components of these 
clusters are shown in[9l Within these results, clusters of A- type monomers are 
found to be either non-percolating, or to be percolating in only one dimension. 
Time series data for the fraction of A-type monomers in clusters of each type 
are shown in[lOl While no obvious energetic relaxation is observed, the CaSPA 
results suggest significant structural relaxation. 

The CaSPA data clearly shows that, across most of the simulation, the 
system is dominated by the presence of a single, large (containing more than 
70% of all A-type monomers) cluster which is percolating in one dimension. 
As such, we can suggest that the system is equihbrating towards a hexagonal 
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Figure 10: Fraction of A-type monomers in non-percolating clusters (/o) and in 
clusters which are percolating in one dimension (/i) as a function of simulation 
time, t. Multiple points at a given time step indicate multiple clusters of that 
type at that time step. 

cylindrical liquid crystal. 

Up until t ~ 19300At, the percolating cluster coexists with one or more 
smaller, though still significant, non-percolating clusters, which appear to be 
micellar in character. Where only a single micelle exists, it can be seen to be 
extended, and approximately cylindrical (Ai w A2 < A3) (see snapshot configu- 
ration in[7l[a)). The largest principle component of this micelle can be seen to be 
of the order of half the length of the simulation box (A3 « 7rc). Where two mi- 
celles exist, the smaller micelle is approximately spherical (Ai « A2 ~ A3 w Tc), 
while the larger micelle is extended and approximately cylindrical (see snap- 
shot configuration in[71^b)). The excess material contained in these micelles is 
absorbed into the larger cluster at t « 19300At. The smallest principal mo- 
ments of gyration Ai and A2 increase at around t w 19300At, corresponding to 
absorption of this material into the large micelle. 

After t « 19300At, the vast majority (greater than 99.5%) of A-type monomers 
exist in a single cluster. Across most of this time (approximately 85%), that 
cluster is percolating in one dimension (see snapshot configuration in [71[d)), 
however, there are five events where that cluster becomes non-percolating (see, 
for e.g., snapshot configuration in[7ljc)). During these events, the cluster keeps 
approximately the same values for the smaller principal moments of gyration, 
while taking on a large value (A3 ~ 18.5rc) for the largest principal moment. 

The character of the large, percolating cluster does not appear quite cyHn- 
drical from the data; the aspect ratio of the cluster is not unity, but instead 
A2/A1 w 1.9. This value does not change significantly with time. The reason for 
this apparent non-cylindrical character is made clear in the "stitched" snapshot 
configurations ([7]); the cluster does not He along a straight path, but is instead 
significantly curved. This suggests that the configuration is not fully equili- 
brated, but is suffering from a finite size effect, where the cylindrical aggregate 
is not aligned properly within the simulation box, and must bend in order to 
maintain self-contact across periodic boundaries. We suggest that the observed 
events where the cluster becomes non-percolating correspond to equilibration 
towards the correct alignment. 
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4 Discussion 



In the previous sections, we have presented CaSPA, a novel algorithm to char- 
acterise the size and orientation of percolating aggregates, and have tested this 
algorithm on a triblock copolymer system. Results using CaSPA show that this 
system forms percolating one-dimensional aggregates. These results demon- 
strate the strengths of the algorithm. Configurational energy does not indicate 
any relaxation phenomena during the simulation ([6]), whereas measurement 
of the size of aggregated clusters clearly identifies structural relaxation 
While the standard Hoshen-Kopelman algorithm would identify the presence 
of this relaxation, CaSPA allows a fuller description by identifying percolating 
aggregates. In the case presented here, this allows very long elongated mi- 
celles to be distinguished from true percolating cylindrical aggregates. Further, 
CaSPA returns the correct principle moments of gyration for branched perco- 
lating structures (such as the percolating aggregates in[7|(a) and (b)). Finally, 
the "stitching" step in CaSPA can be used to generate snapshot configurations 
of aggregates lO which are easier to visualise than the equivalent "unstitched" 
configurations. Coupled with data on the size of aggregates, this clearer visual- 
ization can aid both understanding and communication of observed behaviors 
in simulation of mesostructured materials. 

Use of CaSPA is complementary to other order parameter measurements 
used to describe mesophase structures. While the algorithm provides a broad 
range of information on the size and shape of aggregates, it does not provide 
information on the degree of segregation of aggregating particles from the rest 
of the system. As a comparative example, the P order parameter proposed by 
Groot et al p[2] provides information on the degree of segregation of particles, 
and can differentiate between structures, but is not intended to provide clear ge- 
ometric information on those structures. In summary, CaSPA provides a useful 
addition to methods for identification and characterisation of mesostructured 
materials. 

References 

[1] T. Smart et al, Nanotoday 3, 38 (2008). 

[2] H. Klok and S. Lecommandoux, Adv. Mater. 16, 1217 (2001). 

[3] D. Siegel and R. Epand, Biopys. J. 73, 3089 (1997). 

[4] S. May, D. Harries, and A. Ben-Shaul, Biophys. J. 78, 1681 (2000). 

[5] I. Koltover, T. Sladitt, J. Radler, and C. Safinya, Science 281, 78 (1998). 

[6] V. Ortiz, S. Nielsen, M. Klein, and D. Discher, J. Polym. Sci. B: Pol. Phys. 
44, 1907 (2006). 

[7] M. Gradzielski, Curr. Opin. Colloid In. 9, 256 (2004). 

[8] S. Sakurai et al, Macromolecules 26, 485 (1993). 

[9] R. Groot and T. Madden, J. Chem. Phys. 108, 8713 (1998). 



14 



[10] C. Soto-Figueroa, M. Rodriguez-Hidalgo, J. Martinez-Magadan, and L. 
Vicente, Macromolecules 41, 3297 (2008). 

[11] A. Khokhlov and P. Khalatur, Chem. Phys. Lett. 461, 58 (2008). 

[12] R. Groot, T. Madden, and D. Tildesley, J. Chem. Phys. 110, 9739 (1999). 

[13] J. Hoshen and R. Kopelman, Phys. Rev. B 14, 3438 (1976). 

[14] P. Hoogerbrugge and J. Koelman, Europhys. Lett. 19, 155 (1992). 

[15] J. Koelman and P. Hoogerbrugge, Europhys. Lett. 21, 363 (1993). 

[16] D. Frenkel and B. Smit, Understanding Molecular Simulation, second edi- 
tion ed. (Academic Press, ADDRESS, 2002). 

[17] R. Groot and P. Warren, J. Chem. Phys. 107, 4423 (1997). 

[18] W. Press, S. Teukolsky, W. Vetterling, and B. Flannery, Numerical Recipes 
in C, second edition ed. (Cambridge University Press, ADDRESS, 1992). 



15 



